Investigation of Gene Networks in Three Components of Immune System Provides Novel Insights into Immune Response Mechanisms against Edwardsiella tarda Infection in Paralichthys olivaceus

Simple Summary Paralichthys olivaceus, a common marine fish, is susceptible to various harmful organisms. One such organism, Edwardsiella tarda, is a particularly potent pathogen that can cause severe illness and even death in fish during prolonged infections. The gills, blood, and kidneys, which are crucial components of the immune system in these fish, play vital roles in immune responses, including in the differentiation of immune cells, removal of diseased cells, and other immunity-related processes. In this study, we infected Paralichthys olivaceus with Edwardsiella tarda and analyzed the genetic information of three components at different timepoints after infection. By using advanced analytical techniques, we identified specific genes and gene networks associated with the immune response. Our innovative approach allowed us to gain a deeper understanding of how the immune system of Paralichthys olivaceus responds to Edwardsiella tarda infection in these important components. These findings provide valuable genetic resources for studying immunity in Paralichthys olivaceus and contribute to our knowledge of the molecular mechanisms underlying Edwardsiella tarda infections in fish. Abstract As a quintessential marine teleost, Paralichthys olivaceus demonstrates vulnerability to a range of pathogens. Long-term infection with Edwardsiella tarda significantly inhibits fish growth and even induces death. Gills, blood, and kidneys, pivotal components of the immune system in teleosts, elicit vital regulatory roles in immune response processes including immune cell differentiation, diseased cell clearance, and other immunity-related mechanisms. This study entailed infecting P. olivaceus with E. tarda for 48 h and examining transcriptome data from the three components at 0, 8, and 48 h post-infection employing weighted gene co-expression network analysis (WGCNA) and protein–protein interaction (PPI) network analysis. Network analyses revealed a series of immune response processes after infection and identified multiple key modules and key, core, and hub genes including xpo1, src, tlr13, stat1, and mefv. By innovatively amalgamating WGCNA and PPI network methodologies, our investigation facilitated an in-depth examination of immune response mechanisms within three significant P. olivaceus components post-E. tarda infection. Our results provided valuable genetic resources for understanding immunity in P. olivaceus immune-related components and assisted us in further exploring the molecular mechanisms of E. tarda infection in teleosts.


Introduction
The Japanese flounder (Paralichthys olivaceus) constitutes a notable marine teleost species recognized for its tender meat and high nutritional value. This economically valuable fish is extensively cultivated along the coastal regions of China, Japan, Korea, and Libraries were constructed following the methods described in Li et al. [34]. Sequencing of the libraries was performed using the Illumina HiSeq 4000 platform. Clean reads were obtained by removing low-quality raw reads, reads with >10% unknown nucleotides, and those containing adaptor sequences. TopHat (version 2.0.13) was utilized for mapping clean reads to the reference genome of P. olivaceus. We employed Cufflinks (version 2.2.0) for transcript construction and abundance estimation from TopHat results, quantified as fragments per kilobase of transcript per million mapped reads (FPKM). Using Cuffdiff (version 2.2.0), we statistically compared differentially expressed genes (DEGs) between three groups of samples at the same timepoint in each component. The criteria for differential gene screening were q-value ≤ 0.05 and |log 2 fold change| ≥ 1. We selected the union of DEGs at each timepoint of the three components as the gene dataset for constructing co-expression networks.

Construction of Gene Co-Expression Network
The construction of gene networks was performed using the R package WGCNA following the methodology outlined by Ponomarev et al. [35]. Initially, the DEGs were utilized as input for WGCNA, and the correlation coefficients between each gene pair were calculated. The determination of the optimal power parameter (β) was based on achieving a correlation between genes that adhered to the scale-free topology network, with the criterion of R-squared (R2) exceeding 0.8, as described by Zhang et al. [36]. Subsequently, a hierarchical clustering tree was constructed based on the calculated correlation coefficients. Genes demonstrating similar expression patterns were assigned to the same module while unassigned genes were grouped into the grey module. The thresholds for module clustering were set at a minimum module size of 25 and a cutting height of 0.86.

Identification of Key Modules and Genes
To identify the module most strongly associated with the target trait, we calculated correlation coefficients between modules and traits using module eigengenes. In this study, modules exhibiting the highest positive correlation with the trait were designated as key modules. To elucidate the principal functions of genes within these modules, we performed GO and KEGG enrichment analyses. Genes with high connectivity are more likely to hold significant biological importance. Consequently, we selected the top 20 most highly connected genes from each module as key genes, and among these, the three genes with the highest connectivity were designated as core genes. The interactions between key and core genes were visualized using Cytoscape. Furthermore, we investigated the functional linkages among genes in different modules by constructing protein-protein interaction networks utilizing STRING (version 11.5) [37]. Hub genes, which were identified based on their substantial number of protein interactions, were further validated by their higher connectivity.

qRT-PCR Validation
To evaluate the accuracy of the RNA-Seq results, we conducted quantitative reverse transcription-polymerase chain reactions (qRT-PCRs). The qRT-PCR protocol followed the methodology described by Wang et al. [38].
Total RNA was extracted using the TRIzol method based on the manufacturer's instructions. RNA quality, concentration, and integrity were determined by using 1.2% agarose gel electrophoresis, NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA), and Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA), respectively. RNase-free DNase I was employed to eliminate DNA contamination. A volume of 350 µL of deproteinized solution (RW1) was added to the adsorption column (RA) and centrifuged at 12,000 rpm for 30 s, and the waste liquid was removed. Then, 50 µL DNase I was added to the solution and allowed to stand at room temperature for 15 min, followed by centrifugation at 12,000 rpm for 30 s after adding 350 µL RW1. Then, mRNA was purified using oligo (dT) magnetic beads from RNA and randomly sheared into pieces of approximately 200 base pairs in the fragmentation buffer. First-strand cDNA was synthesized from fragmented mRNAs under the action of reverse transcriptase and random hexamer primers. DNA polymerase I and RNase H was used to synthesize second-strand cDNA. After connecting the adapters and PolyA, PCR was used to isolate the adaptermodified fragments. We validated the expression levels of 21 core or hub genes, and their respective names and primer sequences can be found in Table 1. Gene-specific primers were designed using Primer Premier 5.0 software (PREMIER Biosoft, San Francisco, CA, USA). As a reference gene with stable expression, β-actin was employed. In each group, three biological replicates were used in qRT-PCR. We added 10 ng template cDNA to 20 µL SYBR Premix Ex Taq II (TaKaRa) solution. Then, cDNA was incubated for 5 min, undergoing 45 cycles of 95 • C for 15 s and 60 • C for 45 s in a LightCycler 480. Finally, the melting curve was analyzed, and fluorescent signal accumulation was recorded at the 60 • C, 45 s phase during each cycle by using a LightCycler 480 (Roche, Basel, Switzerland). Additionally, 2 −∆∆Ct method was used for statistical analysis.

Screening of DEGs
Following a comprehensive differential expression analysis, we identified 808 and 1265 DEGs in the blood of P. olivaceus infected with E. tarda at 8 and 48 h, respectively; 456 and 1037 DEGs in the gills at 8 and 48 h; and 1319 and 4439 DEGs in the kidneys at 8 and 48 h. Overall, we uncovered a total of 6459 DEGs (the union of all DEGs in the three components). These DEGs were subsequently employed for further analyses, and their distribution was visualized using a hierarchical clustering heatmap ( Figure 1).

WGCNA
In this study, the 6459 DEGs were employed to construct gene co-expression networks. A scale-free topology model fit was constructed, and a power β value of 6 was selected to ensure that the network held high biological significance ( Figure 2a). We concurrently assessed alterations in the mean connectivity of DEGs at different power β values (Figure 2b). Consequently, we identified twenty modules with module sizes ranging from 20 to 1359 (Table 2, Figure 3).

WGCNA
In this study, the 6459 DEGs were employed to construct gene co-expression networks. A scale-free topology model fit was constructed, and a power β value of 6 was selected to ensure that the network held high biological significance ( Figure 2a). We concurrently assessed alterations in the mean connectivity of DEGs at different power β values ( Figure 2b). Consequently, we identified twenty modules with module sizes ranging from 20 to 1359 (Table 2, Figure 3).

Identification and Functional Analysis of Key Modules
Our investigation led to the identification of six modules exhibiting the strongest correlation with infection-responsive components or infection timepoints. These modules were based on gene eigenvalues. Specifically, the magenta, blue, and green modules demonstrated significant correlation with blood, gills, and kidneys post-infection, respectively ( Figure 4a); DEGs in the yellow, lightgreen, and cyan modules played significant roles in components at 0, 8, and 48 h post-infection (Figure 4b). GO and KEGG enrichment analyses elucidated immune response mechanisms in the three components of P. olivaceus infected with E. tarda. For instance, 77 GO terms (covering biological processes, cellular components, and molecular functions) were enriched in the magenta module while 67, 281, 28, 492, and 149 terms were enriched in the cyan, green, lightgreen, yellow, and blue modules, respectively (Table S1, Figure 5). Several terms, such as protein export from nucleus, NIK/NF-κB signaling, innate immune response, and chemotaxis, were associated with immune responses. Similarly, pathways such as NOD-like receptor signaling, cAMP signaling, chemokine signaling, and the biosynthesis of antibiotics signaling in the six modules substantially contributed to immune response processes (Table S2, Figure 6). 281, 28, 492, and 149 terms were enriched in the cyan, green, lightgreen, yellow, and blue modules, respectively (Table S1, Figure 5). Several terms, such as protein export from nucleus, NIK/NF-κB signaling, innate immune response, and chemotaxis, were associated with immune responses. Similarly, pathways such as NOD-like receptor signaling, cAMP signaling, chemokine signaling, and the biosynthesis of antibiotics signaling in the six modules substantially contributed to immune response processes (Table S2, Figure 6).

Construction of PPI Network and WGCNA Key Networks
We created distinct gene networks for the three modules demonstrating the highest correlation with infection-responsive components, employing 20 key genes from each module (Table S3, Figure 7a). The core genes exhibiting the highest connectivity in the magenta, green, and blue modules were xpo1, herc1, and slc8a1; nop58, ppan, and cpne4; and atp8b1, src, and trim39, respectively. These genes participated in immune responses of different components and played pivotal roles therein. Likewise, three networks correlating with infection timepoints were constructed (Table S4, Figure 7b). Core genes like mefv, rarres2, and loc109626689 were identified as central to the cyan module, while grina, slc12a3, and hcar2 emerged as core members of the lightgreen module. Moreover, ttc36, aqp, and pyroxd2 were identified as core genes in the yellow module. These genes modulated immune responses at three distinct timepoints. Analyses of protein-protein interaction networks across modules revealed that xpo1, src, and dkc1 ( Figure 8a) might play crucial roles in the three infected components, and stat, mefv, and tlr13 ( Figure 8b) could have important effects within the first 48 h of infection. The interaction numbers and connectivity of hub genes are detailed in Table 3.
lated immune responses at three distinct timepoints. Analyses of protein-protein interaction networks across modules revealed that xpo1, src, and dkc1 (Figure 8a) might play crucial roles in the three infected components, and stat, mefv, and tlr13 ( Figure 8b) could have important effects within the first 48 h of infection. The interaction numbers and connectivity of hub genes are detailed in Table 3.   lines represent higher correlation. Interactions between modules are shown by lines connecting genes. (b) Gene networks representing different infection timepoints.

Validation of Hub Genes Using qRT-PCR
To validate the accuracy of our transcriptome analysis results, we performed qRT-PCR experiments. The qRT-PCR data demonstrated concordance with the RNA-Seq findings, with all measured DEGs producing single products (Figure 9). The consistent gene expression trends observed using both RNA-Seq and qRT-PCR methods further support the reliability and accuracy of our RNA-Seq results.

Validation of Hub Genes Using qRT-PCR
To validate the accuracy of our transcriptome analysis results, we performed qRT-PCR experiments. The qRT-PCR data demonstrated concordance with the RNA-Seq findings, with all measured DEGs producing single products (Figure 9). The consistent gene expression trends observed using both RNA-Seq and qRT-PCR methods further support the reliability and accuracy of our RNA-Seq results.

Comprehensive Analysis of WGCNA and PPI Network
The utilization of WGCNA and PPI networks constitutes two cutting-edge meth ologies in the field of transcriptome analysis. These analytical techniques enable the

Comprehensive Analysis of WGCNA and PPI Network
The utilization of WGCNA and PPI networks constitutes two cutting-edge methodologies in the field of transcriptome analysis. These analytical techniques enable the determination of gene interactions and the identification of hub genes most closely associated with the subject matter being investigated. In previous research, these two methods have primarily been employed individually, rather than in conjunction with one another, to identify hub genes [38,39]. In the current study, we not only identified key and core genes with high connectivity within each key module but also innovatively ascertained numerous additional hub genes potentially involved in the regulation of P. olivaceus immunity through a combined approach using PPI network and WGCNA analysis. The results demonstrate that some genes possess both high connectivity and a considerable number of protein interactions, therefore designating them as hub genes likely implicated in the immune response processes in P. olivaceus. Concurrently, even though certain genes exhibit relatively low connectivity, their high number of protein interactions suggests a plausible significance in P. olivaceus immunity. The integrative analysis of these two methodologies allows for the discovery of numerous hub genes imperative for regulating the immune response of P. olivaceus, ultimately providing support for further studies examining the immune response processes in P. olivaceus following E. tarda infection.

Functional Analyses of Key Modules and Genes
We employed WGCNA analysis on 6459 DEGs to identify key genes affecting immune response processes in P. olivaceus, resulting in six modules with 120 key genes that were most relevant to specific components and infection timepoints. GO and KEGG analyses revealed numerous immune-related GO terms and KEGG signaling pathways enriched within various modules, signifying that module-enriched genes were intimately related to immune response processes. Among them, genes in magenta and lightgreen models were only enriched in one signaling pathway each. It was possible that there were too few genes in two modules to be enriched into multiple pathways. It was worth noting that the hub gene xpo1, enriched in the magenta model, was enriched in the only pathway. It involved and regulated biological immune responses processes, suggesting that xpo1 might bear the main immune functions of fish blood. In our examination of the PPI network of 60 key genes in three key modules correlated with infection time and component data, we identified common interactions among these key genes even though they were associated with different infection times and component types. This observation implies that these genes may cooperatively regulate P. olivaceus immunity. By identifying six hub genes with high protein interaction numbers or high connectivity, we were able to posit that these hub genes play significant roles in the immune response processes of the selected components following E. tarda infection. Investigating these hub genes will aid in advancing our understanding of the immune response mechanisms in P. olivaceus after infection.

Analysis of Magenta Module Associated with Blood Immunity
In this module, three genes, namely xpo1, herc1, and slc8a1, were identified as core regulators of immune response processes. Exportin 1 (xpo1), also referred to as chromosomal maintenance 1 (crm1), is a critical transport receptor belonging to the transport protein superfamily. It is ubiquitous in eukaryotic cells and mediates the transportation of numerous proteins and diverse RNA species [40]. Previous investigations have demonstrated that xpo1 can facilitate the transport of several immune proteins, including PI3K, NF-κB, AKT, p53, and Wnt, to various tissues and organs [41,42]. These proteins significantly influence the growth, proliferation, apoptosis, inflammatory responses, and activation of T and B cells, among other processes [43]. Consequently, xpo1 may indirectly modulate immune responses by regulating the transport of immune proteins. Moreover, xpo1 impacts the transportation of growth-regulatory proteins. The upregulation of xpo1 can enhance the transport of growth regulator proteins, thereby affecting immune cell proliferation [44]. In the current study, xpo1 was identified as a core gene due to its high connectivity and interaction with six key genes in modules related to gill and kidney infection. We hypothesize that xpo1 contributes to the growth and proliferation of immune cells within the blood and modulates P. olivaceus gill and kidney immune responses by transporting immune proteins and cells.
Hect domain and RCC1-like domain-containing protein 1 (herc1), a member of the HERC family proteins, regulates various cellular processes including cell cycle progression and cell signal transduction [45]. As a ubiquitin ligase, herc1 primarily controls cell migration through the activation of the MKK and p38 modules [46]. Simultaneously, herc1 modulates the RAF, MEK, and ERK signaling pathways, promoting the growth, proliferation, differentiation, and transport of immune cells [47].
The solute carrier family 8 member A1 (slc8a1) is an intracellular calcium regulator that mediates Na + and Ca 2+ exchange [48]. The upregulation of slc8a1 suppresses the proliferation and migration of pathological and necrotic cells, promoting apoptosis [49]. Previous research has also demonstrated that slc8a1 can regulate biological inflammatory responses [50]. While slc8a1 has been extensively studied in humans, its function in fish remains largely unexplored. The immune roles of slc8a1 in fish blood require additional investigation.
Collectively, these three core genes primarily function as transporters of immune cells and regulators of cell growth, differentiation, and apoptosis. The sustained upregulation of these genes suggests that immune response processes in P. olivaceus infected blood can effectively inhibit pathogen invasion.

Analysis of Green Module Associated with Kidney Immunity
Three core genes primarily related to kidney immunity were found to potentially participate in and regulate various immune response processes. The majority of Box C/D non-coding RNA was located within nucleoli, wherein it interacted with several proteins to create RNA-protein complexes (RNPs). These complexes played a pivotal role in the modification and processing of pre-ribosomal RNA [51]. nop58, a primary interacting partner of Box C/D RNA, facilitated gene methylation upon binding, thus modulating the expression of immune-related genes [52,53].
ppan, a critical biogenic factor in nucleoli and mitochondria, regulated mitochondrial homeostasis [54]. Emerging research indicates that ppan is involved in immune response processes by interacting with the WNT signaling pathway, thus regulating the proliferation and differentiation of immune cells [55]. Furthermore, it contributes to the autophagy process responsible for degrading and recycling infected cells and organelles [56].
These three genes displayed the highest correlation with fish kidney immunity. Among them, nop58 and ppan genes might be integral to immune response processes in P. olivaceus kidneys such as immune cell differentiation and immune-related signaling pathway regulation. Currently, no direct link between cpne4 and immune response has been identified. However, we preliminarily hypothesize that cpne4 may be associated with P. olivaceus kidney immunity. The specific functions of this gene require further exploration in future experiments.

Analysis of Blue Module Associated with Gill Immunity
The gill plays a crucial role in initiating immune responses, specifically in the proliferation and differentiation of immune cells. In this study, we identified atp8b1, src, and trim39 as potential regulators of gill immunity in P. olivaceus.
atp8b1 is a phospholipid transporter that maintains epithelial cell stability via phospholipid transport [59]. Mutations in atp8b1 have been associated with progressive familial intrahepatic cholestasis type 1 (pfic1) in humans and mice [60,61]. However, the function of atp8b1 in fish has not been investigated, and its exact mechanisms remain unclear. We hypothesize that atp8b1 in our study may regulate the production and apoptosis of gill epithelial immune cells through phospholipid transport, but further research is needed to confirm this hypothesis.
src is a non-receptor protein-tyrosine kinase that has been studied extensively in recent years [62]. It is implicated in regulating various cellular processes including proliferation, differentiation, and adhesion [63]. Notably, cell adhesion is closely associated with inflammatory responses and is a key step in immune responses [64]. src also participates in intracellular signaling that modulates the expression of multiple cytokines and immune complexes [65]. Previous research has demonstrated that src regulates the immune responses of immune receptors, C-type lectins, and integrins [65,66]. Salmond et al. [67] found that src influences the activation and differentiation of T cells by regulating the T-cell receptor signaling pathway, thereby modulating immune response processes. Our results indicate that src is enriched in the cell adhesion term, suggesting its significant immune function in the cell adhesion process in fish gills following E. tarda infection. Concurrently, src interacts with five genes enriched in other modules, suggesting its potential role in regulating the proliferation, differentiation, and other immune cellular processes in fish blood and kidneys.
The TRIM protein family plays an essential immunomodulatory role by activating various immune-related signaling pathways [67,68]. As a key immune gene in this family, trim39 regulates the cell cycle and apoptosis of immune cells [68]. Previous studies have demonstrated that trim39 is a significant regulator of innate immune responses in diverse organisms including humans and fish [68,69]. Moreover, trim39 can activate the NF-κB signaling pathway, regulating cellular immunity, inflammatory response, and other crucial immune processes [70].
In conclusion, we posit that the gill, as the first line of defense against pathogen invasion, may play a significant role in P. olivaceus immunity following E. tarda infection. Investigating the functions of these three core genes could enhance our understanding of the immune response mechanisms in P. olivaceus gills.

Analysis of Yellow Module Associated with 0 h Infection
ttc36, aqp8, and pyroxd2 were identified as hub genes enriched in the yellow module, indicative of their potential role in the activation and proliferation of immune cells. ttc36, a member of the tetratricopeptide repeat (TPR) subfamily, has received little attention in previous research [71]. TPR proteins, recognized as scaffolds, have been implicated in numerous biological processes such as cell cycle regulation, apoptosis, protein transport and degradation, and resistance to bacterial infection [72]. In light of these previous findings, one might posit that ttc36 similarly contributes to such functions. Notably, Jiang et al. demonstrated that ttc36 could interact with HSP70 and upregulate its expression, leading to enhanced resistance against pathogen invasion [73]. This suggests a possible role for ttc36 in the regulation of heat shock proteins and the modulation of biological immunity.
Aquaporin 8 (aqp8), a protein responsible for maintaining the intracellular water balance, has recently been implicated in processes such as tumorigenesis, transport, and apoptosis [74]. aqp8 facilitates H 2 O 2 transport and regulates various cellular processes [74,75]. Recent evidence has also revealed a role for aqp8 in B cell signal transduction, wherein it modulates biological immunity through the activation and differentiation of B lymphocytes [76].
pyroxd2, localized within the inner membrane and matrix of mitochondria, is involved in the regulation of mitochondrial function [77,78]. Additionally, pyroxd2 governs cellular processes, selectively promoting apoptosis in diseased cells while sparing their healthy counterparts [77].
These findings collectively suggest that certain P. olivaceus genes exhibit immune functions even in the absence of infection, perhaps as a defense mechanism against bacterial or viral contaminants at relatively low concentrations in seawater. Interestingly, these three core genes predominantly manifest in the kidney, with minimal expression in the blood. This observation implies that the kidney may serve as the primary organ for immune function in uninfected P. olivaceus, a hypothesis requiring substantiation through further experimentation.

Analysis of Lightgreen Module Associated with 8 h Infection
The lightgreen module analysis identified three core genes, namely grina, slc12a3, and hcar2, which are associated with 8 h infection in P. olivaceus. grina, also referred to as trim3, is an anti-apoptotic protein present in cell membranes and has been implicated in the regulation of cell survival and apoptosis, as opposed to the regulation of Ca 2+ concentration and neurotransmitter release [79]. Studies on grina-mediated apoptosis in various organisms such as mice, humans, and zebrafish have been conducted [79,80]. grina has been found to play a crucial role in regulating apoptosis during neuronal development and endoplasmic reticulum (ER) stress [81]. In light of these findings, we speculate that grina might regulate immune cell apoptosis in P. olivaceus following infection.
slc12a3, a sodium chloride cotransporter primarily expressed in the kidney, is regulated by multiple hormones such as aldosterone, glucocorticoids, and insulin [82,83]. The involvement of slc12a3 in complex bio-ion concentration regulation processes implies its role in maintaining the ion balance within an organism [84,85]. Previous research has demonstrated slc12a3's role in regulating Na + concentrations in fish [86], which is further supported by our findings revealing an enrichment of the sodium transmembrane transport signaling pathway.
hcar2 functions as a G-protein-coupled receptor involved in immune and inflammatory responses [34]. By regulating cAMP, NF-κB, and other signaling pathways, hcar2 promotes T-cell activation and modulates cellular processes such as immune cell proliferation and differentiation, thereby influencing immune responses [34,87]. Our study posits that hcar2 might play a similar role in P. olivaceus immunity.
All three core genes displayed a significant upregulation within the 8 h infection period, indicating that P. olivaceus may undergo intricate immune response processes following infection. It is noteworthy that while hcar2's immune functions have not been reported previously, our results highlight a high expression of hcar2 after E. tarda infection. Consequently, we propose that hcar2 might have critical immune functions in fish, warranting further investigation in future studies.

Analysis of Cyan Module Associated with 48 h Infection
In P. olivaceus, intricate immune responses emerged at 48 h post-infection, indicative of the potential involvement of three core genes as key regulatory factors in immune response processes. One such gene, rarres2 (also known as chemerin), functions as an immunomodulatory factor primarily produced in adipose and cutaneous tissues [88]. Chemerin has been implicated in various immune-related processes, including microbial defense and inflammation, functioning as a leukocyte chemokine [89]. For instance, rarres2 can recruit innate immune cells to inflammatory sites, subsequently eliciting inflammation [90]. Additionally, chemerin has been reported to directly modulate cellular processes such as the migration, invasion, and proliferation of pathological cells in humans [91]. Our findings suggest that rarres2 may both be involved in the regulation of immune cell proliferation and differentiation as well as contribute to inflammatory processes.
Another core gene, mefv, has been implicated in regulating inflammatory responses and various cellular processes [92]. The pyrin protein encoded by mefv participates in innate immune processes, producing IL-1β to mediate inflammation [93]. Furthermore, mefv enrichment in the ubiquitin-protein transfer activity signaling pathway implies a potential role in the ubiquitination of immune proteins, thus regulating differentiation and apoptosis [92,93]. Despite extensive research on mefv in humans, its function in fish remains largely unexplored [94]. Consequently, the effects of differential mefv expression in fish warrant further investigation.
In this module, we also identified a core gene, loc109626689, which has yet to be annotated. Given its high connectivity, we postulate that loc109626689 may be closely related to the immune response in fish. Functional analyses of other key and core genes in this module further suggest that loc109626689 could participate in the regulation of inflammatory responses. Interestingly, the core genes exhibited elevated expression levels at 48 h post-infection and performed significant immune functions, excluding the uncertain gene loc109626689. To date, the immune functions of these genes have not been extensively examined in fish. We hypothesize that these genes may exert essential immune-related effects in P. olivaceus such as modulating inflammatory responses and diverse cellular processes.

Gene Function Analyses between Modules
In the present research, a novel and innovative approach was employed, combining PPI network analysis with the conventional WGCNA module connectivity analysis, to jointly pinpoint potential hub genes that could modulate the immune response in P. olivaceus. Our findings identified three crucial hub genes, namely xpo1, dkc1, and src, which appear to regulate immunity in P. olivaceus. dkc1 encodes the dyskerin ribonucleoprotein, implicated in promoting telomerase synthesis and maintaining telomere length and integrity, which in turn propels cell proliferation [95]. Concomitantly, it has been shown to regulate the proliferation, differentiation, senescence, and death of telomere-related cells [96]. Recent research has revealed the role of dkc1 in controlling the proliferation, migration, and invasion of tumor cells [97]. Based on these findings, we tentatively posit that dkc1 may regulate cellular processes such as the proliferation and migration of immune cells in P. olivaceus. xpo1 and src were identified due to their higher protein interaction numbers and higher intramodular connectivity. Following our analysis, it is reasonable to suggest that these genes may exert a significant influence on the immune response processes in P. olivaceus. stat1, tlr13, and mefv were pinpointed as regulators of the immune response processes in P. olivaceus, specifically within the first 48 h post-infection, based on interaction analyses. stat1, a key transcription factor affecting the JAK-STAT signaling pathway, modulates an array of cellular functions including proliferation, differentiation, migration, and apoptosis [70,98]. It has been implicated in various immune processes such as activating T cells, NK cells, and other immune cells and regulating innate immunity, adaptive immunity, and certain inflammatory responses [99]. Concurrently, it fosters the proliferation and activation of microglia to regulate neuroimmunity [100]. Toll-like receptors (TLRs) serve as pattern recognition receptors (PRRs) that are prevalent in a wide range of immune cells, participating in the regulation of immune response processes [101]. They have the capacity to recognize a broad spectrum of pathogens, encompassing bacteria, viruses, and parasites [102]. tlr13, a crucial member of the TLR family, plays an indispensable role in biological immune functionality. It has been shown to activate NF-κB, AP-1, type 1 interferon, and other immunoregulatory factors, thereby modulating innate immune responses [103]. In addition, it promotes the proliferation and differentiation of T cells, B cells, NK cells, neutrophils, and other immune cells, thus influencing biological immune response processes [101,104]. Remarkably, recent studies have demonstrated that tlr13 is expressed in fish kidneys, potentially playing a pivotal role in fish kidney immunity [105]. Based on these conclusions, it is evident that these two hub genes perform essential functions in the immune response processes, leading us to speculate that they may regulate the proliferation and differentiation of immune cells, inflammatory response, the activation of immune factors, and other mechanisms in P. olivaceus infected with E. tarda. Moreover, mefv, another hub gene with higher protein interaction numbers, is speculated to partake in the inflammatory response to modulate immune response processes in P. olivaceus.
Our study indicates that the three component-associated hub genes were enriched in three distinct modules, signifying that the three components of P. olivaceus' immune system may exhibit comparable sensitivity to E. tarda and could all perform vital immune functions to resist and eliminate bacteria following infection. An additional three hub genes were identified in modules correlated with 8 and 48 h infection, potentially acting as hubs to participate in and regulate the immunity of P. olivaceus after E. tarda infection. The immune functions of these six genes within P. olivaceus remain ambiguous; therefore, further exploration of their roles will be valuable for understanding and characterizing the immune response processes in P. olivaceus.

Conclusions
In this study, we employed an innovative approach combining WGCNA and PPI network analysis to investigate the immune response mechanisms of P. olivaceus infected with E. tarda. Through comprehensive analyses, we identified six key functional modules. Within these modules, we identified nine core genes (xpo1, herc1, slc8a1, nop58, ppan, cpne4, atp8b1, src, and trim39) that were closely related to the immunity of the three examined components. Genes regulated unique immune response mechanisms in different components, indicating that fish component immunity may have a high specificity. Additionally, we discovered mefv, rarres2, loc109626689, grina, slc12a3, hcar2, ttc36, aqp, and pyroxd2 as regulators of P. olivaceus immunity at different timepoints. The results indicated that E. tarda may inhibit ion transport, disrupt the ion balance, and promote cell apoptosis in P. olivaceus in the early stages of infection and induce severe inflammatory response during long-term infection.
Furthermore, our PPI network analysis, constructed by key genes in each group, identified six hub genes (xpo1, dkc1, src, stat1, tlr13, and mefv) that displayed high protein interaction numbers and connectivity. These hub genes are likely to have crucial immune functions in P. olivaceus when infected with E. tarda, regulating the proliferation and differentiation of immune cells and inducing inflammatory reactions to resist E. tarda infections. Notably, this research represents the first attempt to integrate PPI network analysis with WGCNA in order to identify multiple hub genes that potentially regulate the immunity of P. olivaceus.
The findings of our study lay the foundation for a deeper understanding of the immune response mechanisms in teleosts and hold promise for addressing the challenges of reduced production and disease in P. olivaceus under high-density artificial breeding conditions. By providing insights into the intricate interplay between genes and immune processes, our research contributes to the advancement of knowledge in the field of teleost immunity and has implications for the development of targeted strategies to enhance disease resistance in aquaculture settings.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani13152542/s1, Table S1: GO terms in each module; Table S2: KEGG signaling pathways in each module; Table S3: Connectivity within modules and protein interaction numbers between modules in Figure 7a; Table S4: Connectivity within modules and protein interaction numbers between modules in Figure 7b.

Conflicts of Interest:
The authors declare no conflict of interest.